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NUMERICAL STUDY OF HYDROGEN-AIR SUPERSONIC COMBUSTION 
BY USING ELLIPTIC AND PARABOLIZED EQUATIONS 


By 

Tawit Chitsomboon 1 and Surendra N. Tiwari 2 
ABSTRACT 

The two-dimensional Navier- Stokes and species continuity equations are 
used to investigate supersonic chemically reacting flow problems which are 
related to scramjet-engine configurations. A global two-step finite-rate 
chemistry model is employed to represent the hydrogen-air combustion in the 
flow. An algebraic turbulent model is adopted for turbulent flow calcu- 
lations. The explicit unsplit MacCormack finite-difference algorithm is 
used to develop a computer program suitable for a vector processing 
computer. The computer program developed is then used to integrate the 
system of the governing equations in time until convergence is attained. 

The chemistry source terms in the species continuity equations are 
evaluated implicitly to alleviate stiffness associated with fast chemical 
reactions. The block mono-diagonal system of algebraic equations, resulting 
from treating the source terms implicitly, can be inverted efficiently on 
vector processing computers. Both prefixed and nonpremixed reacting flows 
are solved. Some of the results are compared with those using a complete- 
reaction chemistry model. It is found that the finite-rate model predicts 
the same trends as the complete-reaction model, but the former allows lesser 
extent of combustion. 

The problems solved by the elliptic code are re-investigated by using a 
set of two-dimensional parabolized Navier-Stokes and parabolized species 

Graduate Research Assistant, Department of Mechanical Engineering and 
Mechanics, Old Dominion University, Norfolk, Virginia 23508. 

2 Eminent Professor, Department of Mechanical Engineering and Mechanics, Old 
Dominion University, Norfolk, Virginia 23508. 

ii 



equations. The thermo-chemical and turbulence models utilized are the sane 
as used in the fully elliptic procedure (with proper modifications). A 
linearized fully-coupled ful ly-impl icit finite-difference algorithm is used 
to develop a second computer code which solves the governing equations by 
marching in space rather than time, resulting in a considerable saving in 
computer resources. Results obtained by using the parabolized formulation 
are compared with the results obtained by using the ful ly-ell iptic equa- 
tions. The comparisons indicate fairly good agreement of the results of the 
two formulations. 


i i i 


FOREWORD 


This is a progress report on the research project "Analysis and Compu- 
tation of Internal Flow Field in a Scramjet Engine." The period of perform- 
ance on this research was December 1, 1985 through May 31, 1986. The work 
was supported by the NASA Langley Research Center (Computational Method 
Branch of the High-Speed Aerodynamics Division) through research grant NAG- 
1-423. The grant was monitored by Dr. Ajay Kumar of the High-Speed Aero- 
dynamics Branch. 
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1. INTRODUCTION 


General background information and motivation for the present study are 
discussed in this section. Review of related literature is presented fol- 
lowed by specific objectives of the study. 

1.1 Background 

A comprehensive research program has been underway particularly at the 

NASA Langley Research Center to develop supersonic combustion ramjet (scram- 

★ 

ject) engine for hypersonic aircraft cruising at Mach 5 and beyond [1-4], 
Flow characteristics at these speeds dictate engine-airframe integration 
where the forebody of the vehicle is used to partially compress incoming air 
for the inlet of the engine which is mounted directly underneath and becomes 
an integral part of the vehicle. The aftbody of the vehicle, immediately 
downstream of the engine, is used as part of the nozzle to provide addition- 
al thrust. Hydrogen has been a strong candidate for fueling this type of 
engine due to its many desirable properties such as high specific impulse 
and cooling effect. 

The engine module has a rectangular capture area with sweep of the 
sidewalls and a cut back cowl in order to provide flow spillage which allows 
the inlet to start at low flight speeds. The compression process in the 
inlet is completed by wedge-shaped struts which also provide multiple planes 

for injection of gaseous hydrogen fuel. The combustor has a diverging shape 
and employs varying amounts of parallel and perpendicular fuel injection 
from the sidewalls and the struts. 

A cut-and-try procedure which is inevitable in the initial phase of 
experimentation has always proved to be expensive and time consuming 


*The numbers in brackets indicate references. 




especially for hypersonic flow research; this has led the way to numerical 
experimentation. In this study, the attention is directed on applying nu- 
merical techniques to compute a chemically reacting flow field through the 
scramjet engine. 

Numerical investigation of a complex chemically reacting flow requires 
enormous computer resources, well-proven numerical algorithms and realistic 
chemical kinetics models. The advent of supercomputers such as the VPS-32 
vector processing computer at the NASA Langley Research Center has allevi- 
ated some restrictions on the computer resource, provided the computer code 
is written to take advantage of the hardware of the computer. Advancements 
in nunerical algorithms in the past decade have made numerical investi- 
gations more reliable research tools than ever before. Chenical kinetics of 
the hydrogen-air system are the most studied and understood phenomena, pro- 
viding ample reliable data base for all levels of sophistication. 

The flow field of a scramjet engine is quite complex. Shock waves 
emanating from the sidewalls and the struts leading edges are strong enough 
to separate the turbulent boundary layers in the inlet. Interactions of 
shock waves and expansion fans are also expected to be strong. The flow 
field near transverse fuel injectors is particularly complex. Blockage and 
deflection of the supersonic mainstream flow by the fuel jets give rise to 
an adverse pressure gradient which separates the turbulent boundary layer 
upstream of the injectors. The flow field in the combustor is further com- 
plicated by the highly turbulent mixing and exothermic reaction of hydrogen 
and air. Moreover, the efffect of combustion in the combustor could be felt 
upstream through the subsonic boundary layers which could degrade the per- 
formance of the inlet. Several injector-diameters downstream of the injec- 
tor, however, the flow becomes less complex; a predominant flow direction 
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can be determined in which there is no flow separation. 

To resolve all kinds of interaction discussed, the fully elliptic form 
of the governing partial differential equations must be used in the inlet 
and in the near field of the injectors. Flow in the rear section of the 
combustor and aftbody of the vehicle (where flow additionally expands) can 
be solved by a much less expensive parabolic or parabolized form of the 
governing equations since elliptic effects in the predominant flow direction 
are negligible in these regions. 

1.2 Survey of Pertinent Literature 

Nonreacting flows in various two- and three-dimensional inlet configu- 
rations were studied by Kumar [5-7], In these studies, the fully elliptic 
form of the governing equations are employed for accurate representation of 
the flow field which involves separations .due to shock-wave/boundary-layer 
interaction. The computer codes developed are highly efficient on the 
vector processing computers since they are fully vectorized. Limited 
success was reported by Chitsomboon and Tiwari [8] and Chitsomboon et al . 

[9] in using a parabolized form of the governing equations for solving non- 
reacting flow through the related configurations. Effect of gaseous fuel 
injection without combustion and with a complete combustion model were also 
investigated [10-11]. 

Incorporation of chemical reaction into a numerical scheme for solving 
fluid flow problems is quite an endeavor. First of all, a chemistry model 
compatible with the physical problem being solved must be selected. Depend- 
ing on the rate of chemical reaction relative to the characteristic time of 

the fluid flow, the suitable chemistry model could be any one of the follow- 
ing models: 
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1 . 


frozen flow model 


2. finite rate model 

3. equilibrium model 

4. complete reaction model 

In general, the finite-rate model is the most accurate one but it is also 
the most difficult to implement. More often than not, a compromise must be 
made as to the level of complexity of the chemistry model selected and the 
computer resources available. 

Simple introductions to the physics and the mathematics of chemically 
reacting flow can be found in [12-15]. Over the past decade, a number of 
finite rate chemistry models for the hydrogen-air system have appeared in 
the literature. Rogers and Schexnayder [16] proposed as many as 60 reaction 
paths in their model; this is certainly one of the most complete representa- 
tions of hydrogen-air reaction. Unfortunately, the enormous number of reac- 
tion paths and chemical species involved in the model do not lend themselves 
feasible for a numerical investigation, at least within the current capa- 
bility of state-of-the-art computers. Intermediate- level models [17] reduce 
the reaction paths down to twenty-five and eight with the number of species 
of twelve and eight, respectively. Except for some inaccuracies during the 
ignition delay period, the eight-reaction model performs as well as the 25- 
path model. Although these models are considerable less tedious than the 
60-path model, they are expected to be too costly to use for routine para- 
metric studies especially if the fully elliptic governing equations were to 
be used. The global two-step chemistry model of Rogers and Chinitz [18] 
offers an alternative to nunerical investigators since it is a finite rate 
model with only four species of reactants and products and with only two 
reaction paths. This model was deduced by fitting the temperature history 
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of a 28-reaction model [16] used in a series of constant-pressure stream- 
tube calculations. While there are a number of limitations to this model 
such as ignition phase inaccuracy and a tendency to overpredict the flame 
temperature, the model is considered to be appropriate for initial para- 
metric study of overall mixing and extent of combustion. 

It is the rule rather than the exception that incorporation of a finite 
rate chemistry model into a numerical scheme will yield a stiff systan of 
governing differential equations. Stiffness arises from disparity in time 
scales of the governing equations and has its meaning only in the context of 
numerical or computational mathematics. From a mathematical view point, a 
system of equations is stiff if the eigenvalues of the coefficient matrix of 
its locally linear representation have negative real parts that are widely 
disparate [13], From a practical perspective, the size of the time step of 
integration is severely limited by stability rather than accuracy when stiff 
phenomenon exists [19]. The disparity in time scales, in turn, is caused by 
wide differences in rates of chemical reaction of various species; usually 
tne fastest reaction is the one that causes stiffness. 

During transient period (or ignition delay period), some free radicals 
may be produced very rapidly for a very short duration of time; thereafter, 
they stay relatively unchanged at their equilibrium values since other re- 
actions begin to consume them as soon as they are produced. It is paradoxi- 
cal that stiffness occurs only when the rapid transient state has subsided 
[13, 20], During the rapid transient period, small time steps must be used 
anyway in order to resolve the rapidly-changing temporal history and there- 
fore, the systen is not stiff since time step is limited by accuracy. Dur- 
ing equilibrium period, rates of reactions slow down by several orders of 
magnitude and it is natural to use much larger time step (while still 
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maintaining accuracy) for economical reasons. Use of large time steps, 
however, will result in numerical instability because stability of the 
governing differential equations' is still dictated by the smallest time 
scale. An eigenvalue analysis of a model problem clearly illustrates this 
point (see, e.g. , Ref. 21) . 

Stiff phenomenon is studied extensively in the framework of ordinary 
differential equations (ODE) which are representative of the governing equa- 
tions of a well-stirred reactor or a batch-chemistry system (systems with no 
spatial gradients) [19, 22-26], In most cases, some form of implicitness is 
used to evaluate chemistry source terms in order to alleviate stiffness. 

Some of the most popular stiff ODE solvers are those of Hindmash [23-24]. 
These solvers are based on the methods of Gear [25-26] which, in turn, con- 
tain different implementations of the Adams formulae. When applied to flow 
problems, an Adams-based algorithm was found to be inferior to the four-step 
Runge-Kutta algorithm [27]. 

In light of the above discussion, it is realized that fully explicit 
numerical methods cannot be used to time integrate the stiff governing equa- 
tions of chemically reacting flow due to the prohibitively small time step 
required to ensure stability. If only a steady state solution is sought, 
implicit methods could be used with no restriction on time-step size. Use 
of fully implicit methods, however, requre large computer memory and the 
inversion of the block multidiagonal system of algebraic equations. These 
are difficult to implement to take advantage of the hardware of the vector- 
processing computers. An alternative is to evaluate the chemistry source 
terms in the species equations implicitly while other terms, not contri- 
buting to stiffness, are evaluated explicitly. This semi-implicit method 
was proposed by several investigators [28-31]. 
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A stability analysis [31] reveals that by linearizing the implicit 
source terms one ends up with a preconditioning matrix which modifies the 
unsteady terms of the governing equations. The preconditioning matrix 
allows each species to evolve at its own characteristic time scale. With 
implicit source terms, the governing parameter is now the fluid dynamics 
time scale which, in turn, is dictated by the Courant-Friedrichs-Lewy (CFL) 
condition. 

The partially implicit method and the global two-step chemistry model 
[18] was used successfully by Drummond, et al . [32] to calculate premixed 
H 2 -air reacting flow using a spectral method. Bussing and Murman [31, 33 j 
used a finite volume method to solve for flows over ramps and rearward fac- 
ing steps. Chitsomboon and Tiwari [34] used a finite difference method for 
similar problems. Eklund et al . [27], studied relative merits of an Adams 
scheme and a Runge-Kutta scheme. A nonpremixed calculation of this partic- 
ular type of reacting flow was done by Chitsomboon et al. [35]; in the 
study, gaseous hydrogen fuel was injected perpendicular to the sidewalls and 
the fuel struts of an inlet-combustor configuration. The results reported 
in these studies have indicated that the semi-impl icit method together with 
the global two-step chemistry model provides a viable tool in numerical 
study of H 2 -air reacting flow by using the elliptic governing equations. 

As mentioned earlier, the flow field far downstream of the fuel injec- 
tors is relatively less complex; a predominant flow direction exists in 
which there is no flow separation. If only a steady state solution is 
sought, a class of approximation to the fully elliptic equations, para- 
bolized approximation, is physically reasonable. The motivation in doing 
this kind of approximation is that an N-dimensional problem reduces to an 
(N-l)-dimensional problem. Furthermore, the reduced dimension problem can 
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be solved by marching in the flow direction once (or, at most, a couple of 
times) instead of marching in time for many steps as in the elliptic proce- 
dure. This method is highly efficient with regard to memory requirement and 
total CPU time of the computer program. A brief discussion of this partic- 
ular approximation is given below. 

The boundary layer approximation introduced by Prandtl at the turn of 
the century has enjoyed a wealth of analyses both analytically and numeri- 
cally (see, e.g.. Refs. 36-38 and references therein). An order of magni- 
tude analysis shows that, for a high Reynolds number flow, the streamwise 
diffusion terms and the whole normal momentum equation can be dropped since 
they are terms of high order. The remaining steady state equations, with 
the streamwise pressure gradient specified, are mathematically parabolic and 
can be solved by integration in the streamwise direction provided appropri- 
ate initial and boundary conditions are supplied. In its original incep- 
tion, the boundary layer equations were used to correct for thin viscous 
effects near solid walls of the outer inviscid solution. The streamwise 
pressure gradient can be obtained from the inviscid solution; this method 
has been termed the direct method. Due to the neglect of the normal momen- 
tum equation, the direct method cannot resolve the viscous-inviscid inter- 
action which could be important for some classes of flow. To partially 
account for the interaction, modified methods such as the inverse method or 
the interacting boundary layer method must be used. Numerical algorithms of 
these methods could be quite tedious since they require adaptive adjustments 
of various parameters. 

Parabolized approximation contains all the boundary layer equations 
plus the normal momentum equation. The inclusion of the normal momentum 
equation makes the equations very versatile. First of all, the entire Euler 
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equations are included; this eradicates the need of laborious patching of 
viscous and inviscid solutions required by the boundary layer procedure. 

The viscous-inviscid interactions in the normal direction are naturally 
represented and the equations still enjoy the efficient procedure of space 
marching in obtaining the solution. Like the boundary layer approximation, 
this approximation is valid for high Reynolds number flows only. Tradition- 
ally, the equations are termed "parabolic" if the streamwise pressure gradi- 
ent is prescribed a priori. The equations are called "parabolized" if 
mathematical treatments are done to the pressure gradient term such that the 
term can be retained as part of the solution without causing any ill posed- 
ness. Some investigators prefer to call both formulations as reduced 
Navier-Stokes equations (RNS). 

It appears that Ferri [39 j and Rudman and Rubin [40 ] are among the 
first to propose using this class of approximation to the full Navier-Stokes 
equations. Since then, it has found numerous applications from the incom- 
pressible regime to the hypersonic regime for both external and internal 
flow. 

For imcompressible and subsonic flows, characteristics of downstream 
flow can affect those at upstream since the pressure waves which propagate 
disturbances in the flow travel faster than the local fluid velocity. For 
supersonic flow, information can still propagate upstream through subsonic 
boundary layers close to solid walls. In essence, a disturbance at a point 
in the domain affects the solution through out the domain; this constitutes 
an elliptic boundary value problem. If the untreated parabolized equations 
are marched with the pressure determined as part of the solution, then, the 
equations are illposed since they represent an elliptic boundary value 
problem with a parabolic initial value problem. The mathematical 
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characteristic of the equations is controlled by the streamwise pressure 
gradient term in the streamwise momentum equation. Lighthill [41] demon- 
strated the existence of an exponentially growing solution (called the 
departure solution) which reflects illposedness by retaining the streamwise 
pressure gradient term. Lubard and Helliwell [42] showed that for a back- 
ward difference on the streamwise pressure gradient term, the departure 
solution can be avoided if the integration step size (ax) is greater than 
(ax) where Ax is of the order of magnitude of the subsonic layer thickness. 
This is a rather paradoxical conclusion since one would expect that the 
smaller the integration step size the more accurate and the more stable the 
solution will be. However, if ax is less than Ax (which is the extent of 
upstream interaction) the parabolized equations will try to represent 
upstream influence and will result in a departure solution. If the pressure 
field can be prescribed a priori, however, the equations are mathematically 
parabolic and there is no difficulty in marching in the streamwise 
direction. The initial prescribed pressure field can be updated once the 
entire domain is swept. This is accomplished by a Poisson equation for 
pressure derived from the momentum equations. Schemes based on this proce- 
dure are quite expensive since the Poisson equation might have to be solved 
several times [43-45], but it is still relatively economical as compared to 
using the fully elliptic equations. 

Algorithms which are valid from the incompressible regime to the hyper- 
sonic regime which do not require solution to the Poisson equation have been 
reported in [46-48]. These algorithms used forward difference on the pres- 
sure gradient term which necessitates a known or guessed pressure field. 

The departure behavior has been claimed to be taken care of by the forward 
difference procedure. Also, global iterations on the pressure field, by 
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multiple sweeping of the domain, become necessary in order to correct the 
guessed pressure field. Another widely used algorithm was developed by 
Patankar and Spalding [49]. This is a single sweep method for incompress- 
ible flow. The assumed pressure field and the velocity field are corrected 
simultaneously once an integration step is completed. The correction is 
made in a manner such that the continuity equation, which is not included in 
the implicit solution procedure, is satisfied. An extension of the algo- 
rithm to compressible flows was made in [50]. The modified algorithm does 
not seen to capture shock wave well due probably to the nature of the pres- 
sure correction procedure [51]. Global iteration on the pressure field was 
also studied [52]; it was found that the iteration had quite an effect on 
the characteristic of the flow field. 

For supersonic high Reynolds number flow without streamwise separation, 
the algorithm of Vigneron et al. [53] is very efficient and accurate. The 
accuracy is partially attained by solving all the governing equations in a 
fully coupled manner. The algorithm is really a steady state version of 
those proposed by Beam and Warming [54 J and Briley and McDonald [55]. Since 
it is a non-iterative single sweep method, the algorithm is very efficient. 
The only cumbersome part of the method is that one has to evaluate the 
Jacobian matrices of all the flux vectors and solve them in a block tri- 
diagonal manner instead of just a scalar tri-diagonal procedure encountered 
in the previously discussed methods. The departure behavior of the solution 
is suppressed by retaining only a fraction, a>, of the streamwise pressure 
gradient within the subsonic portions of a boundary layer. The parameter, 
a), was obtained from an eigenvalue analysis of the locally linearized gov- 
erning equations such that the equations become well posed. By just manipu- 
lating the governing equations, Khosla and Lai [56] have also revealed the 
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existence of the parameter a>. A similar fully coupled fully implicit algo- 
rithm was developed also by Schiff and Steger [57]. An elaborate eigenvalue 
analysis revealed again that the i 1 1-posedness of the system of equations is 
caused by the streamwise pressure gradient in the subsonic layer. To render 
the equations well posed, a sublayer approximation was then proposed. In 
this approximation, the pressure gradient within the subsonic zone close to 
a solid wall is evaluated at the first supersonic mesh point away from the 
wall. In an independent study, the sublayer model was also proposed by 
Rubin and Lin [58]. 

Another good feature of this type of algorithm is that it can be cast 
into a fully conservative form, thus, capable of capturing strong shock 
waves. Rakich [59] and Chitsomboon et al . [9] applied global iteration 
procedures on the pressure field to this algorithm for attached and sepa- 
rated flows, respectively. A modified form of the FLARE approximation [60] 
was used in solving flows with streamwise separations. 

Parabolic governing equations are quite popular in solving supersonic 
chemically reacting flows of hydrogen and air [50, 52, 61-62], The algo- 
rithm used in these studies is that of Partankar and Spalding [49] with a 
compressibility correction for the pressure-velocity update procedure. 
Equilibrium chemistry models were used in all of these studies. Finite-rate 
chemistry models were employed by Rogers and Chinitz [18] for solving simi- 
lar flows in which gaseous hydrogen is injected parallel to the flow direc- 
tion. Some difficulties were reported in using a finite-rate chemistry 
model. Some of the species concentrations became negative and the computer 
program could not continue but there was no discussion as to how the chemis- 
try source terms were evaluated [18]. A finite rate chemistry model was 
also utilized by Sinha and Dash [63]. In their study, the species equations 
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were decoupled from the fluid dynamics equations and the chemistry source 
terms for the next integration step were obtained by a shooting technique. 

As discussed earlier, the Patankar-Spalding-based algorithm probably could 
not capture strong shock waves typical of flow in a scramjet; thus limiting 
the application of the algorithm to relatively simple flows. For flows with 
embedded strong shock waves, particularly those with strong interactions 
between chemistry and fluid, an alternative approach is warranted. 

1.3 Objectives 

The objectives of this study are to, first, incorporate the finite-rate 
chemistry package of Rogers and Chinitz [18] into an existing nonreacting 
elliptic computer code. This code was developed initially for solving flows 
in inlets of scranjet engine configurations [5, 64J. The unsplit MacCormack 
finite difference algorithm was used to advance the discretized form of the 
two-dimensional Navier-Stokes equations in time until a steady state solu- 
tion is attained. The code is fully vectorized and is operational only on 
the VPS-32 vector processing computer at the NASA-LaRC (an upgraded version 
of the CDC Cyber-205). The incorporation of the chemistry package has 
included four species continuity equations along with the original set of 
the governing equations. The energy equation is also modified to make it 
appropriate for a chemically reacting flow. Stiffness associated with fast 
chemical reactions is alleviated by evaluating the chemistry source terms 
implicitly. The resulting computer code will be used to solve various re- 
acting H 2 -air problems related to scramjet engine configurations. 

A second computer code, based on fully implicit parabolized metnods of 
Vigneron et al . [53] and Schiff and Steger [57], is also developed. The 
reason for using this particular algorithm is to capture strong shock waves 
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typical of the flow problems for which the code is developed. The thermo- 
chemical models employed here are the same as those used in the fully 
elliptic procedure. The chemistry equations and the fluid dynamics equa- 
tions are solved in a fully implicit fully coupled manner in order to 
resolve the fluid-chemistry interactions. The results obtained by using 
this space marching code will be compared with those obtained by using the 
more expensive elliptic code. 

The basic formulations for this study are described in Section 2. The 
eliptic and the parabolized governing equations are presented together with 
the thermo-chemical models and other important formulations. The finite 
difference algorithms for both the elliptic and the parabolized equations 
are presented in Section 3 where other related informations are also pro- 
vided. The important results obtained in this study are presented and dis- 
cussed in Section 4. Finally, specific conclusions and some suggestions for 
further study are provided in Section 5. 

2. BASIC FORMULATION 

This section deals with the discussion of the physical model and basic 
governing equations for the elliptic and parabolized formulations. The 
relations for the chemistry and thermodynamic models are also provided. 

2.1 Physical Model 

As pointed out in the introduction, the primary goal of this study is 
to compute chemically reacting flow through the scram jet engine. A 
schematic of a scramjet engine is shown in Fig. 2.1 in which it is shown 
that the engine consists of many identical modules; the size of a module is 
appropriate for ground testing. An engine module and a part of its topview 
are also illustrated. The fuel is injected from the fuel struts as 
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indicated by the arrows. Figure 2.2 depicts the complex flow field in the 
vicinity of a fuel injector. To resolve the complexity of the flow in the 
near field of the injector, the fully elliptic form of the governing 
equations must be used. The parabolized form of the governing equations as 
discussed in the introduction can be used in the region far downstream of 
the fuel injector in which the flow becomes less complex. 


2.2 Elliptic Governing Equations 

The two-dimensional Navier- Stokes equations and the species continuity 
equations in body-fitted coordinate systan written in the strong 
conversation law form can be expressed symbolically as 


+ 3E + 3F = W 
3t 3? 3n J 

! where 


( 2 . 1 ) 


q = i [p, pu, pv, pH-p, p, ] T 
J 


( 2 . 2 ) 


E = 


pu 


PUU + Vxx - Vxy 
pvu + y n T yx - x n T yy 

('>H-PW n (UT x /VT xy *q x )-x n (VT yy 4-UT yx *q y ) 

P,u + y„(B( f ,) x > - x„(S<f 1 ) y ) 


(2.3) 
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Fig. %*2 Flow field near an injector 



F = 


pv 


puv - y T + XT 

5 xx 5 xy 

pvu - y_x w + x _ t 

5 yx 5 yy 


(BH-p)v-y 5 (uT xx tVT xy *q x )tx 5 (uT y)< *VT yy +q y ) 

p^v - y c (e(f i ) x ) + x 5 (e(f 1 ) y ) 


(2.4) 


Here, x^ denotes 3x/3£, and so forth, and 


u = y u - xv 
n n 

(2.5a) 

v = - y^u + x ? v 

(2.5b) 

J * U 5 y n - v*)" 1 ' 

(2.5c) 


It is to be understood that the source terms, W, for the Navier-Stokes 

equations are identically zero. For the species equations in source terms 

represent rates of species production or extinction to be discussed in 

detail in the subsequent section. The quantities t , t , and t are 

xx- xy yy 

components of the stress tensor. With the Stokes hypothesis assumed, these 
quantities are given by 


_ „ / . w 4 3u 2 3 Vn 

r xx - P - (» + uJ ( ) 

1 33x 33y 


• ( U * W t ) (— * 11) 

3y 3x 
_ _ . , , x (2 3u 4 3v x 

3 3x 3 3y 


T = T 

xy yx 


yy 


(2.6a) 

(2.6b) 

(2.6c) 


The quantities q and 


q^ are components of the heat 


flux; the fluxes are due 
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to heat conduction and transport of enthalpy by species diffusion. For a 
laminar flow the terms are given by 


= - k 11 - 

N 

P 01 

3f . 
[(— 

ax 

1-1 

ax 

= - k£ 

N 

- pD E 

af . 

[(— : 


1=1 

ay 


where 


hi = no + / T Cp . di . 


(2.7a) 


(2.7b) 


( 2 . 8 ) 


Note that the effective binary diffusion coefficient, D, is used for all 
species. Together with the assumption that the Lewis number is unity, the 
approximation provides a great simplification to these heat fluxes. The 
assumption, while not always stated explicitly, appears to be a fairly rou- 
tine assumption in the analysis of gas phase reacting systems. With the 
manipulation presented in Appendix A, the heat flux terms can be written 
simply as 


* 


V 


_ _ j 3h 


Pr 3x 


- C— ) 


ah 


Pr ay 


For a turbulent flow the terms are modified as 


(2.9a) 

(2.9b) 
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Q x = - (— + — ) — (2.10a) 

Pr Pr t 3x 

q y = - (— + — ) — • (2.10b) 

J Pr Pr^ 3y 

The molecular viscosity, u, is assumed to be dependent upon temperature and 
is calculated by the Sutherland's formula: 

* I 1 ' 458 * 10 ~ 6 > T3 “ ■ (2.11) 

T + 110.4 

Pure air properties are used in the above relation since the mixture is 
dominated by the nitrogen species. The turbulence eddy viscosity, u^, is 
obtained by the algebraic two-layer eddy viscosity model of Baldwin and 
Lomax [65]. In this model, the inner layer is represented by the Prandtl- 
Van Driest mixing length formulation and the Clauser approximation is used 
in the outer layer. The model is very convenient to use since there is no 
need to find the edge of a boundary layer like some other algebraic models. 
The study of Peters et al. [66] shows that it is a viable turbulence model 
and provides a good starting point in turbulence modeling. The laminar and 
turbulent Prandtl nimbers are assumed to be equal to constant values of 0.72 
and 0.9, respectively. 

The subscript i's in the species equations range from one to four and 
denote the species 0 2 , H 2 0, H 2 and OH, respectively. The mass fraction of 
nigroten can be obtained by the mass constraint relation: 
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( 2 . 12 ) 


f = i - Z f . 
b 1-1 1 


The quantity 3 in the species diffusion terms is defined as 


a « - (-!!— + — ) 

Le*Pr Le t *Pr t 


where 


(2.13) 


Le ■ a/D » 1 


(2.14) 


and 


Pr = v/o . 


(2.15) 


The laminar and turbulent Lewis mmbers are assimed to be equal to one. 

2.3 Parabolized Governing Equations 
The two-dimensional parabolized Navier-Stokes and parabolized species 
equations in a body-fitted coordinate system are expressed in a strong 
conservation law form as 


where 



ff + ^ F 1 + F V> 
8? 3n 


W. 

i 

J 


(2.16) 


- [0, E (l-a»)p , E (l-ai)p, 0, 0, 0, 0, 0] 


(2.17) 
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puv 

puv 

5 y 

+ __ 

pVV + up 

pHu 

j 

pHv 

. P i U - 


V . 


pu 


pv 

puu + p 


PUV 

puv 

♦ 'j. 

pvv + p 

pHu 

J 

pHv 

. P i U . 


. P i V - 


F v ■ 


XX 


xy 


ut + vt + q 
xx xy H x 

(6f i } x 


1 

J 


0 

T 


y* 


yy 

ut + vt + q 

yx yy 


(6 Vy 


The derivatives in the viscous flux, F, are given as 


T 


XX 


(u + u t ) 




\> 


T xy " - * “t> <\ \ + \ v n> 

T yy ' (l,+ “t ) ( ^ n y v n ' \ \ u n> 


(2.18) 

(2.19) 

. ( 2 . 20 ) 

(2.21a) 

(2.21b) 

(2.21c) 
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(2.22a) 


q = - ( JL + _L) n x h 

x Pr Pr , x n 


q y = - (— + — ) n h 
" Pr Pr t ^ n 


(2.22b) 


(6 f i) x = n (Bf .) 
i x x in 


(Bf.) = n (Bf.) 

ry y rn • 


(2.23a) 

(2.23b) 


The above parabolized equations are obtained by dropping the unsteady terms 
and the- streamwise viscous derivative terms. For high Reynolds number 
flows, neglect of these terms should not degrade quality of the solution by 
much since they are high order terms. Note also that only a fraction u of 
the streamwise pressure gradients are retained. The parameter w is a quan- 
tity of magnitude varying between zero and one and is obtained from an 
eigenvalue analysis [53]. The relation for a> is given in terms of the 
streamwise Mach nunber as 


u = YM x 

1+(y-1)MJ ' 


(2.24) 


From the above relation, it can be seen that <o is equal to zero when the 
streamwise Mach number, M , equals zero; it is equal to unity when M equals 

A A 

one. For M x greater than one, i.e., supersonic flow, u> is assigned the 
value of one. The parameter w is necessary in order to suppress upstream 
pressure interactions. From a mathematical viewpoint, to suppresses the 


23 


ellipticity of the equations rendering them well posed for marching in the 5 
direction. It is natural to expect that this technique will work well only 
for supersonic high Reynolds nunber flows where subsonic layers close to 
solid walls are thin. For flows with large subsonic zones, alternative 
methods, such as those discussed in the introduction, are recommended. 

In the supersonic region, the full Euler equations are included; the 
equations should have capability to capture strong shock waves. The normal 
pressure gradient is also retained in full; this allows free pressure inter- 
action in the cross stream direction without having to resort to the explic- 
it interaction mechanism used in bound ary- layer techniques. 

The quantities y, Pr and 8 are obtained in the same manner as those for 
the elliptic equations. The vorticity which is required in the algebraic 
turbulence model is evaluated by retaining only the cross stream component. 


2.4 Chemistry Model 

A chemistry model is needed in order to obtain the chemistry source 
terms, W, as functions of the dependent variables, q. In this study, the 
global two-step finite-rate chemistry model of Rogers and Chinitz £18 j is 
adopted. The model is given by 


k 


H2 + O2 


2 OH 



(2.25) 


2 OH + H 2 



(2.26) 


where the k 's are the forward reaction rate constants and the k ,'s are the 
f b 
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backward reaction rate constants. The chemistry source terms can be obtain- 
ed by applying the law of mass action to the two reaction equations. For a 
general reaction equation: 


f 

J T i 

Z A. . C. 


j=l 


ij J 


E B. . C. j i 1,2,. ..,1 
j=l 1J J 


(2.27) 


the law of mass action states that the rate of change of molar concentration 
of the (j)th species by the (i)th reaction is given by 


(C'.) . = (B. . 
J i ij 


-V 



J 

tt 

k=l 




(2.28) 


The total rate of change in molar concentration of the ( j ) t h species is 
obtained by sunming over all the reaction equations 


I 

c;. = z (c;). (2.29) 

J i=l J 1 


The source terms on mass basis, W, are finally found by multiplying the 
molar changes by the corresponding molecular weights 


W.=C;M.. (2.30) 

J J J 

By applying the law of mass action to the present two-step global model, the 
chanistry source terms on the mass basis are given as 
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k fl P 1 p 3 


(2.31a) 


W 2 = - 


k bl M 1 p I 


w 2 = 


2 k f2 M 2 p 3 p| 


2 k 


b2 2 


M 3 1< 


m 2 


W 3 = M 3 W! . m 3 w 2 

Mi 2M 2 


W 4 = - 2 M4 Wi _ M^ 
Mi M 2 


where subscripts 1,2,3, and 4 are used to represent 0 2 , H 2 0, H 2 , and 
respectively. 

The forward reaction rates are computed from the Arrhenius' law 


f i 


N. -E./R°T 
= A i T e 


For the present model the rates are given by 


K fl = A x 


T -10 e -4865/R°T 


k 


f2 



T -13 e -42500/R°T 


where 


Ai = (8.9170 + 31.433/0 - 28.95) (10) 44 m 3 /kmol-sec. 


(2.31b) 

(2.31c) 
(2.3 Id ) 
OH, 

(2.32a) 

(2.32b) 

(2.32c) 

(2.33a) 
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A 2 = (2 + 1.333/0 - 0.8330) (10) 58 


m 6 /kmol 2 -sec. 


(2.33b) 


R° = 1.1972 cal/mole-K. 


(2.34) 


Note that the pre-exponential constants, A.., are dependent upon equivalence 
ratio, the functions were determined by curve fitting such that the rates 

agree well with the more sophisticated model [18]. The backward rate con- 
stants are related to the forward rate constants by 


k. . * k /K 
di f i i 


(2.35) 


where K^'s are equilibrium constants. The values for K/s are given by 


Kj = 26.16 e 


■8992/T 


(2.36a) 


69415/T 

K 2 = (2.682 x 10-9)(T) e m 3/kmol . (2.36b) 


It should be noted here that the global two-step chemistry model is not 
accurate for predicting flames with long ignition delay times for tanpera- 
ture on the order of 1000 K. Specifically, the reaction would proceed con- 
tinuously without a delay period even if the mixture temperature is below 
the ignition temperature. To prevent this nonphysical phenomenon, a cut off 
temperature of 1000 K was used to control the reaction, i.e., if tanperature 
of the mixture is less than 1000 K, reaction is not permitted. 

The model is used unmodified for the parabolized formulation. 
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2.5 Thermodynamics Model 

In general, thermodynamic properties of the mixture are computed by 
sunming properties of individual species weighted by species mass fractions. 
The specific heat at constant pressure of the ( i ) th species is assumed to be 
a linear function of temperature 


C . 
pi 


a. + b.T 
i i 


(2.37) 


where a., and b.. are constants which are obtained by curve fitting the therm- 
ochemical data of Ref. 67. The numerical values of these constants are 

given in Table 2.1. The static enthalpy of the mixture can be expressed as 

I 

i 


5 

h = E f.H. + C T (2.38) 

i=l 11 p 

where 


5 

« E f. (a- + 0.5b. T) . (2.39) 

P i=1 i i i 

The total enthalpy can now be given by 

H = h + 0.5 (u 2 + v 2 ) . (2.40) 

The mixture gas constant can be obtained by a mass weighted summation over 
all species as 


5 

T? = E f.R. . (2.41) 

1=1 
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Table 2.1 Numerical values of various constants 


Species 

H 8 (Joule/kg) 

a 

b 

R( Joule/kg-K) 

0 2 

-271,267 

0.1198 

947 

254 

h 2 o 

-13,972,530 

0.4391 

1858 

457 

h 2 

-4,200,188 

2.0546 

12867 

4035 

OH 

+1,772,591 

0.1656 

1673 

478 

n 2 

-309,483 

0.1035 

1048 

240 


29 



The equation of state for the mixture can be written as 


p * pRT (2.42) 

Finally, the specific heat ratio of the mixture is calculated from 

C 

Y = (2.43) 

V 1 

3. METHODS OF SOLUTION 

In this section the finite difference algorithms for the elliptic and 
the parabolized formulations are presented. Also discussed are the boundary 
conditions and the smoothing functions used in both formulations. The grid 
generation technique is given in section 3.5. Some important procedures in 
obtaining solutions are finally described in section 3.6. 

3.1 Algorithm for Elliptic Equations 
The symbolic form of the elliptic equations in body-fitted coordinate 
is written here again as 


*■-(!*♦*- Q) (3 

3t 3C 3n 

where Q = W/J for the species equations. Equation (3.1) can be forward 
differenced in time as 
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(3.2) 


q P+1 


q n = - At (ii" + 1!"- q n+1 

35 9n 


where (n+1) denotes the time level where solution is being sought. Note 
that the source terms, Q, are evaluated at the implicit time level (n+1) to 
alleviate stiffness associated with fast chemical reactions as discussed 
earlier. Since the source terms are nonlinear functions of the dependent 
variable vector, they must be linearized to render linear discrete equations 
for solution. By using the Newton linearization scheme, one obtains 

Q n+1 > Q" ♦ ( 12 .)" (q" +1 . ,") . (3.3) 

3q 


The dependent variables, q, in the above equation are the p^/'J of the four 
species. Since both Q and q are vectors of four elements, the term 3Q/3q 
becomes a 4x4 Jacobian matrix. Upon defining: 


M = 3Q/3q 

n n+1 n 
Aq = q - q 

R = 1£ + — - Q 
35 3n 

and substituting Eqs. (3. 3-3. 6) into Eq. (3.2) one obtains 


(3.4) 

(3.5) 

(3.6) 


( I —A t M) n Aq n = - At(R n ) 


(3.7) 
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where I is the identity matrix. Equation (3.7) yields explicit relations 
for the Navier-Stokes part whereas the species continuity equations give 
rise to a block mono-diagonal system of algebraic equations. Components of 
the Jacobian matrix, M, are given in Appendix B. The unsplit MacCormack 
algorithm [68] can now be applied to Eq. (3.7) as follows: 

Predictor step: 


(I -At M) n Aq n = - At r£ 

(3.8a) 

„n+l _ „n , . n 
q = q + Aq 

(3.8b) 

Corrector step: 


( I -At M) n+1 Aq n = - At r£ +1 

(3.9a) 

q n+1 = q n + 0.5( Aq n + 

Aq n ). (3.9b) 


The subscripts f and b in Eqs. (3.8a) and (3.9a) denote forward and backward 
finite difference, respectively. 

3.2 Algorithm for Parabolized Equations 
For convenience, the symbolic form of the parabolized equations in 
body-fitted coordinates are expressed here again as 

!£+ ?I + — (F t + FJ = Q. (3.10) 

n 3n 
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The Beam and Warming finite difference algorithm [54] is used to discretize 
Eq. (3.10). The algorithm was initially proposed for time dependent equa- 
tions; the modifications for the present steady state equations are obvious 
upon replacing the marching direction with C coordinate. The algorithm for 
Eq. (3.10) can then be written as 

n_ L (AC)a r 3 n . n„-, 

A E + 2 — — [_ A (F t + F v ) - A Q] 

1 + b 3n 

= - A? [— (Fj + F ) - Q] n + — - — A n 1 E 
3n 1 v 1 + b 

- — +0(a-b - 0.5) ( AC ) 2 + 0(AC) 3 - (3.11) 

35 

Many numerical schemes can be generated by varying the two parameters, a and 
b. Some familiar schemes are given below: 
a = 0 , b = 0 ; Euler explicit 

a = 1 , b = 0 ; Euler implicit 

a = 1 , b = 0.5 ; Three-point backward 

a = 0.5, b = 0 ; Trapezoidal 

If the difference between a and b is equal to one half, the algorithm is 
second order accurate in the C direction, otherwise it is first order accu- 
rate. 

Since the flux terms and the chemistry source terms are all nonlinear 
functions of the dependent variable vector and the metric quantities of 
coordinate transformation, they need to be linearized. Unless an appropri- 
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ate linearization scheme is adopted, the conservative property of the algo- 
rithm given in Eq. (3.11) will be disrupted. Conservative linearization is 
very important particularly for internal flow where mass flow rate across 
the channel has to remain constant. Also, shock waves will not be captured 
accurately if the linearization procedure disturbs the conservative property 
of the algorithm. Another important consideration is that the linearization 
scheme should honor the formal accuracy of the algorithm. In this study, a 
scheme very similar to that used by Schiff and Steger [57] is used to lin- 
earize the nonlinear terms. The general technique is to express the terms 
to be linearized as functions of the dependent variable vector, q, and the 
metric quantities; only the terms involving q are linearized, the metrics 
are frozen at the implicit level. This is nothing but a basic mathematical 
operation for a function of multiple variables. The nonlinear terms are 
linearized as follows: 


A n E = E n+1 - E n = [E n + (— )Vq] - [ e"" 1 + + 0(A£) 3 (3.12) 

3q 9q 


rKj 

4 n F, v * (— I>v)Vq + 0(45) 2 (3.13) 

’ 3q 

A n q = (i£)Vq + 0(ac) 2 (3.14) 

3q 


where a tilde over a quantity means that quantity is to be evaluated by 
using q at the indicated level whereas the metric quantities are evaluated 
at the indicated level plus one. Note that the (n+1) and the (n) terms of 
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the streamwise flux, E, are linearized in the same fashion to ensure the 
conservative property in the £ direction; it is also second order accurate. 
Even if a first order accurate algorithm in the £ direction is desired one 
must still linearize in this manner to ensure a conservative scheme. The 
cross stream flux, however, is linearized nonconservatively in £ direction 
since it is already represented in a divergence form in the n direction. 

The first order accuracy in linearization of F and Q is made second order 
accurate since the terms have A£ in multiplication. The chemistry source 
terms are also linearized nonconservatively since they are non-conservative 
source terms in the first place. For convenience in algebraic manipula- 
tions, the following definitions are adopted: 


3E 

= a 

9q 


(3.15) 


— ,v = B t (3.16) 

3q 


^ = C. (3.17) 

3q 


Substituting Eqs. (3.12-3.17) into the original algorithm, Eq. (3.11), and 
rearranging one obtains 
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(3.18) 


[A + L J1 f_i (8 »B J - C))Vq = - a"- 1 ? 

1+b 3n 

+ A n ~V _1 q - [J. (Fj+F v ) - Q] n + — A n-1 E 

1+b 3n 1 v 1+b 

- — + 0 (a-b-0.5) (A?) 2 + 0(ac)3. 

The quantity in the brackets on the left-hand side of the above algorithm is 
an operator that applies on Aq n rather than a multiplication per se. It is 
appropriate to note here that the dependent variable vector, q, in the pres- 
ent formulation is not the same one as used in the elliptic formulation. 

The intrinsic variable for the energy equation is selected as T/J rather 
than (pH-p)/j. The choice of this variable is quite arbitrary. The selec- 
tion was made because it permits easy evaluation of the Jacobian matrices; 
it also gives a good diagonal dominance property for the momentum and the 
energy equations. As a matter of fact, many choices of dependent variaole 
vectors were studied. They included: 


q s [p. 

u, v, T, p.] T 

(3.19) 

q = [p. 

pu, pv, pH, p 1 ] t /j 

(3.20) 

q s [p, 

pu, pv, pT, p 1 .] T /J 

(3.21) 

q = [p. 

PU, PV, T, P^/J • 

(3.22) 


It was found that the last choice, Eq. (3.22), gave the best results. For 
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nonreacting flow, it has been traditional to use 


q = [p, pu, pv, pe t ] T /J (3.23) 

as the dependent variable vector [9, 53, 57]. The use of these variables, 
however, is not convenient for reacting flows since the pressure is related 
to other variables in a complicated way. 

The linearized algorithm given in Eq. (3.18) yields a block tri-diago- 
nal system of algebraic equations; the blocks are 8x8 matrices. By select- 
ing the order of the dependent variables as in Eq. (3.22), the diagonal term 
of the Jacobian matrix. A, for the continuity equation is always equal to 
zero. Experiences with nonreacting flow studied in the past [8, 9] have 
indicated that the algorithm is not very robust; some kinds of implicit 
smoothing functions and/or stabilizing functions are always needed to pre- 
vent a division by zero during the inversion of the block tri-diagonal 
matrices. In this study, a method is proposed to improve robustness of the 
algorithm. In this method, the streamwise velocity in the continuity equa- 
tion is partially lagged one step behind the implicit level where solution 
is being sought. To illustrate the point, the streamwise flux in delta form 
neglecting the metric quantities is given as 

A n E = (pu) n+ ^ - (pu) n - (3.24) 

The flux is then partially lagged as follows: 

A n E = C c u n p n+1 + (1-C c )(pu) n+1 - C c u n “V - (1-C c )(pu) n (3.25) 
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where C c is a parameter of magnitude less than one. This modification 
renders the diagonal term of the continuity equation to be C c u n rather than 
zero. Although the species continuity equations do not yield zero diagonal 
terms, their streamwise velocity is partially lagged in the same way as the 
continuity equation to ensure compatibility. For a second-order accurate 
scheme, the lagging procedure causes the continuity and the species 
equations to be somewhere between first and second order accurate depending 
on the magnitude of C c . The conservative property, however, is not 
disturbed since the terms at (n+1) and at (n) are lagged in the same 
manner. 

Some difficulties arise in evaluating the Jacobian matrix of the 
viscous flux. A typical component of the viscous flux is given as 

— — f, (3.26) 

3n 8n J 3n J 



This can be manipulated further to give a familiar-looking central differ 

i 

ence scheme for a second order derivative. 

i 
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If the flow to be investigated is of high Reynolds number and largely 
supersonic, a single sweep method will, in general, yield a very good solu- 
tion. In this method, one has to march through the domain in the 5 direc- 
tion only once. The term 3P/3? can be included by lagging it one step 
behind the step where the solution is being sought; this could result in a 
departure solution since the elliptic effect is introduced for all practical 
purposes. For a departure-free solution the term 3P/3£ should be dropped 
[53]. The accuracy loss in neglecting this term should not be severe due to 
the thinness of the subsonic boundary layer. Reasonable results for non- 
reacting internal flow, using this method, were reported by Chitsomboon and 
Tiwari [8J. For flows with large subsonic pockets and/or streamwise sepa- 
rations, however, a multi-sweep technique becomes necessary. The first 
sweep is used to establish a pressure field for the entire domain. The 

term 3 P/3? is then used in the subsequent sweeps to relax and update the 
pressure field until a convergence criterion is met. Usually, the criterion 
is selected by the requirement that the maximum change in pressure from one 
sweep to the next be less than a specified tolerance for all the mesh 
points. The method is also known as the global pressure relaxation tech- 
nique. When flow reversal in the streamwise direction is encountered, a 
simple procedure is to use the FLARE approximation of Reyhner and Flugge- 
Lotz [60]. A negative streamwise velocity (i.e., separated flow) will 
render the marching procedure to become unstable since the governing equa- 
tions are again ill posed. The eigenvalue analysis in Ref. 57 indicates 
that when a flow reversal occurs the positive viscosity tends to amplify a 
disturbance instead of damping it. In the FLARE approximation one simply 
replaces the negative convection velocity with zero or a small positive 
number. The procedure works well only for flows with small separation 
bubbles; even if this is so, an additional procedure is required to make the 
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multi-sweep method stable [9]. Some results of the multi-sweep technique 
for nonreacting internal flows are shown in Ref. 9. For the present 
reacting-flow investigation, only the single sweep method will be used. 

It should be noted that Eq. (3.18) is appropriate for a two-dimensional 
problem. For a three-dimensional problem, an approximate factorization 
technique may be needed in order to avoid inverting a sparse block matrix 
system of algebraic equations. 

3.3 Boundary and Initial Conditions 
Boundary and initial conditions are required for both elliptic and 
parabolized equations; they are discussed in the following subsections. 

3.3.1 Conditions for Elliptic Equations 

Before starting the time integration, initial distribution of dependent 
variables for all interior mesh points are specified using free stream 
values. This choice is very convenient since it requires no computation 
whatsoever and seems to work well for all problems investigated in this 
study. All characteristic values of the supersonic governing equations are 
positive indicating that the inflow conditions should be fixed whereas free 
outflow conditions could be used. In light of this, the inflow conditions 
are fixed at constant free stream values at all time steps. A method of 
extrapolation is employed at the outflow boundary. In this method, depend- 
ent variables at the exit boundary are equated to the variables at one 
station upstream of the exit, for a zeroth order extrapolation; first and 
second order extrapolations can be devised also. The method appears to work 
well even though it violates the rule of the method of characteristics with- 
in the subsonic layer close to a solid wall; this is probably because of the 
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tninness of the subsonic layer. At a solid wall, the no-slip conditions (u 
= v = 0) along with the adiabatic wall condition (3T/3y = 0) or specified 
wall temperature condition are imposed. Based on the analysis of boundary 
layer equations, the zero normal pressure gradient (3p/3y = 0) is also used 
at the solid wall; close examinations of numerical output reveals that this 
is indeed a good approximation. The non-catalytic conditions are used at 
the solid wall for all four of the species continuity equations. In this 
method, normal gradients of species mass fractions are set equal to zero 
since the condition requires no net diffusion of species into the wall. In 
the case where gaseous hydrogen fuel is injected from solid surfaces, the 
conditions of injected hydrogen are specified at mesh points within which 
lie the injectors. At a plane of symmetry, reflection boundary conditions 
are imposed. 

3.3.2 Conditions for Parabolized Equations 

Unlike the elliptic procedure, initial conditions (at the inflow bound- 
ary) have to be as accurate as possible. Any inaccuracy in initial condi- 
tions will jeopardize the accuracy of the solution through out the domain 
since the equations are steady and the method of solution is direct, without 
any relaxation. Any method can be used to specify the initial conditions 
insofar as it provides accurate variable profiles. In this study, the 
initial profiles are obtained from the results of the elliptic code. This 
may seem too restrictive since an elliptic problem has to be solved first 
before a parabolized investigation can be made. In real applications, the 
elliptic equations will be used only in the near field of the fuel inject- 
ors; the exit conditions of the elliptic results are then used as starting 
profiles for the parabolized code. 
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At a solid wall, the conditions imposed are the same as those used in 
the elliptic equations discussed in the previous subsection except that they 
have to be expressed in implicit delta forms compatible with the mmerical 
algorithm. This is another advantage in using T as the intrinsic variable 
for the energy equation; if this is not so, the linearization of variable 
may be necessary in obtaining the implicit boundary conditions. Outflow 
conditions, however, are not needed since the solution is obtained by for- 
ward marching in the flow direction. 

3.4 Numerical Dissipation Functions 
It is well known that the central finite difference schemes, as used in 
this study, exhibits spurious numerical oscillations across a shock front. 
Without any control, the oscillation can destabilize the numerical scheme. 

To suppress this high frequency oscillation phenomenon, fourth-order dissi- 
pation functions are added to the righthand side of the elliptic algorithm 
(Eqs. 3. 8-3. 9) as follows: 


f 4 = UO 4 — — (u + a) [Cj fif P + C 2 6|T] — 
5 9£ J * * H 


(3.28) 


where 


f 4 = Un) 4 l-JL (v + a) [C!«2 p + C 2 S2T] ±- q 


(3.29) 


6?P = 


P i+l,j ~ 2P i,J + P i-l,j 
P i+l,j + 2P i,j + P i-l,j 


6?T = 


T i,j+1 ~ 2T i,J + P i,J-l 
T i ,j+l + 2T i,j + 


(3.30) 


(3.31) 
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The form of the functions were suggested by MacCormack [69]. In its origi- 
nal form, only the term 6 2 P was used; Drummond [70] found that for some 
problems this is not sufficient and suggested inclusion of the term <5 2 T as 
well. The coefficients Ci and C 2 must be selected by numerical experiment. 
For the cases investigated in this study, both coefficients are fixed at a 
constant value of one-half; this value appears to work equally well for 
other problems. Finite difference representations of these smoothing terms 
follow the predictor-corrector sequence of the numerical algorithm. 

For the value of the parameter C equals to zero, the parabolized code 
cannot be marched even for a single station since the diagonal term of the 
continuity equation is identically zero. To render the term non-zero as 
well as to provide some damping effects, the following third order implicit 
function 


-Cj (An) 2 (V n A n )A n q 

is added to the lefthand side of Eq. (3.18). A fourth order explicit func- 
tion given by 

-C A (An) 4 (v A ) 2 q n 
t n n 

is also added to the righthand side of the algorithm. The explicit term 
plays a major role in damping the oscillation. A simplified Von Neumann 
analysis shows that both Cj and must be positive (but note the minus 
signs in front of the coefficients). The value of must not be larger 


A3 


than 1/8 for stability reasons. For large value of Cj, however, the stabil- 
ity bound of increases to about one-half of C^. An algebraic relation 
for which optimal damping of high frequency oscillation is attained is given 
by 


1 + 4 Cj - 16 C e = 0 (3.32) 

It is suggested that be selected first and determined such that the 
above relation is satisfied. By numerical experiments it is found that = 
1/4 and = 1/8 gives the best results at least for the cases considered in 
this study. Different smoothing funtions for this type of algorithm can be 
found in [9, 71, 72]. It should be noted that the explicit function has the. 
Jacobian matrix. A, in multiplication whereas the implicit function does 
not. The Jacobian matrix is included to mimic the form of the linearized 
streamwise flux. Inclusion of the Jacobian matrix in the implicit function 
would destroy its diagonal-dominance effect for which the term is added. 
Numerical experiment indicated that Jacobian scaling of the explicit func- 
tion does give a better smoothing effect than without the scaling. 

3.5 Grid Generation 

An algebraic grid generation technique similar to that developed by 
Roberts [73] is used to obtain computational grids. By using the grid gen- 
eration technique, a physical domain is transformed into a rectangular com- 
putational domain with uniform mesh spacing. The size of the computational 
mesh is arbitrary; in this study, it was selected to be unity. Working with 
the rectangular computational domain is much simpler than with the skewed 
physical domain since finite difference representations of derivatives and 


44 


impl enentation of boundary conditions become fairly straight forward without 
having to resort to interpolations. 

The technique allows clustering of mesh points near boundaries in order 
to resolve high gradient quantities. Clustering can be made either at one 
end or both ends of the coordinate. An algebraic stretching function can be 
given as: 

S . M + < N2 ~ [(6+2°) ° - 8+2»] (3.33) 

2 (2«+l) (1+®) 

where N1 is the smallest coordinate index and N2 is the largest coordinate 
index. If the parameter a is equal to zero the mesh near N1 will be clus- 
tered; the mesh near both N1 and N2 are clustered if a is equal to one half. 
The degree of clustering is controlled by the magnitude of 8. The closer 
the magnitude of 8 is to unity from above the more the clustering will be. 
The quantity a is a variable defined by 


a 


(0 + 1 ) 
(8 - 1 ) 


e-g 

1-a 


(3.34) 


where 0 is a variable ranging from zero to one and is given as 


N-Nl 

N2-N1 


(3.35) 


where N ranges from Nl to N2. 

The magnitude of the physical coordinate at each mesh point can now be 
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obtained by the relation (the normal coordinate, y, is given here for illus- 
trative purpose): 


Y N = [Y N2 (S N ' Nl) + V N2 ' V ^ /(N2 ' N1) (3>36) 

where the subscript N denotes the index of the mesh point in question. The 
stretching of the physical coordinate in the streamwise direction, if de- 
sired, can be achieved in a similar fashion. No attempt has been made to 
enforce orthogonality of the grid system. Orthogonality is not very criti- 
cal since the geometries of interest do not vary appreciably in the normal 
direction. 

The mesh spacing near a solid wall is governed by the Reynolds number 
of the flow. In order to properly resolve the details within a laminar 
boundary layer, the minimum mesh spacing close to the solid wall should 
approximately satisfy the following relation [69] 


AY 


min 



(3.37a) 


For a turbulent boundary layer the relation is given by 


AY min 


100 L 


Re 


0.9 


(3.37b) 


This is achieved by varying the parameter B and can, at best, be a trial and 
error process. 

Figure 3.1 shows a 31x31 grid system for a problan solved in this 
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3.1 


Grid structure of case No. 1. 



study. It can be seen that only the normal coordinate is clustered. Clus- 
tering of mesh point in both normal and streamwise directions is illustrated 
in Fig. 3.2. This is a 87x60 grid used in the other problem that is inves- 
tigated. 


3.6 Solution Procedures 

Some important procedures used in the computer codes for obtaining the 
solution are discussed in this section. 

3.6.1 Procedures for Elliptic Code 

The elliptic computer code is written specifically for the VPS-32 
vector processing computer at the NASA Langley Research Center. The fluid 
dynamic equations are integrated in time from an initial solution by the 
fully explicit unsplit MacCormack finite difference algorithm. The species 
equations give rise to a block mono-diagonal system of algebraic equations 
which can be solved very efficiently on the VPS-32 computer. The LU decom- 
position is applied on these blocks one element at a time. Elements of 
matrices which have the same species indices are traced through spatial 
indices vectorially throughout the domain. Without this technique, the code 
was quite expensive to use (about one hour of CPU time on the VPS-32 comput- 
er). With the present technique implemented, however, the total CPU time 
was reduced to about 800 seconds on the same test case whereas mixing along 
(without combustion) required some 600 seconds. The vectorized subroutines 
that solve the block diagonal system are given in Appendix C. 

Once all the conservative variables are obtained at a time step, the 
primitive variables have to be extracted. Some difficulties arise in deci- 
phering the temperature, T, since it is embedded in pH-p, the intrinsic 
dependent variable of the energy equation. This study assumes that the 



specific heat, and the gas constant, IT, of the mixture can be lagged one 
time step without causing any significant error. With the above assumption, 
temperature can be obtained as follows 

T = [(pH-p) - Ip(u 2 + v 2 ) - e f.H’J/tpr* - PR*) (3.38) 

2 i=l p 

where the starred quantities are to be evaluated at the old time level. 8y 
using Eq. (3.38) the temperature of the new time level can be obtained 
directly without having to solve a second order algebraic equation in T. 

The time step used is determined from the CFL condition: 

At = (ct) min [-^ — , — ] (3.39) 

|u|+a |vj+a 

where a is the CFL number of magnitude less than one. Note again that use 
of this large time step is possible only if the chemistry source terms are 
evaluated implicitly. If an accurate temporal history is desired, one must 
choose At in a different manner (see, e.g.. Ref. 21). Experiences gained 
during the numerical experimentation have suggested that the value of o 
should be small in the early phase of time integration. The magnitude of a 
is then gradually brought up to 0.9 in about 100 time steps. Attempts to 
use o = 0.9 at the very beginning yielded negative temperature which was 
probably due to a relaxation process that was too abrupt. 

Unlike the conventional central difference scheme, the present numer- 
ical scheme lacks a decisive, so called, residual to be monitored for a 
convergence criterion. A steady state solution is determined by monitoring 
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the following criteria: 

1. The maximum change in density for all mesh points from a time step 
to the next be less than a specified tolerance. 

2. The L -norm of the change in density from a time step to the next 
be less than a specified tolerance. 

3. The total integration time is equal to or greater than the time 
required for a fluid particle with free stream velocity to travel three 
times the length of the physical domain. 

The tolerances for the first two criteria are usually selected as numbers of 
three to four orders of magnitude smaller than the changes of corresponding 
quantities over the first integration step. Most importantly, the numerical 
solution must be examined in detail by appropriate measures for physical 
plausabil ity. 

3.6.2 Procedures for Parabolized Code 

The parabolized code was developed in a standard F0RTRAN4 language. 
Initial solution profiles are needed to start the program. The profiles 
were obtained from solution of the elliptic code as discussed earlier. For 
the present two-dimensional problems, the storage requirement of the program 
is equivalent to that of a one-dimensional problem. Ironically, the coding 
chore is more laborious than the elliptic code couterpart due to its fully- 
implicit fully-coupled procedure. Two levels of data are needed to obtain a 
solution at the next station. Only one level of data is used to start up 
the program; this is done by using a first order schene by specifying a=l 
and b=0. Attempt has been made to start up with two levels of data by a 
second order scheme; however, this results in instability due probably to 
the fact that the two level initial data contained too much elliptic infor- 
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mat ion. 


Unlike the elliptic code, the temperature, T, is computed directly; 
this eradicates the need to decipher T from a related variable. Some dif- 
ficulties arise in evaluating the Jacobian matrix of the chemistry source 
terms since the kinetic rate constants are strong and complicated functions 
of temperature. The difficulties are alleviated by lagging the temperature 
of the kinetic rate terms; this does not seem to affect the quality of the 
solution. 


The block (8x8) tri-diagonal system of algebraic equations are solved 
by a direct non-iterative procedure. The so-called Thomas algorithm in 
conjunction with an LU decomposition procedure is employed in inverting 
these block matrices. There is no convergence criteria for the code since 
it is a single sweep method. If the situation warranted, a multi-sweep 
procedure could be implemented very easily. It is recommended that the 
nunber of sweeps be limited to only a couple of times; too many sweeps will 
certainly make the code less attractive. Fortunately, the physical problems 


for which the present parabolized fon 


mu i a l IUII 


is applicable (supersonic high 


Reynolds number flows) almost always guarantee that only a couple of sweeps 
are indeed necessary [9], 


4. RESULTS AND DISCUSSION 

Results of the elliptic code are presented first followed by the 
results of the parabolized code. 

4.1 Results of Elliptic Code 

To gain some confidence with the elliptic code, a simple geometry 
problem was solved first. The flow conditions and tne geometry of tne test 
case are depicted in Fig. 4.1. A 31x31 grid system as shown in Fig. 3.1 was 
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used in this test case. This is a premixed case with an equivalence ratio 
(<j>) of unity. The inflow temperature of the mixture was arbitrarily select- 
ed as 900 K in order to study the effect of the explicit ignition switch 
incorporated in the program. As mentioned earlier, the reaction is not 
allowed to take place when the temperature is lower than 1000 K; it is to be 
expected, then, that reactions should occur only behind the shock wave and 
within the adiabatic boundary layers where temperature goes beyond 1000 K. 
Figures 4.2 and 4.3 show the concentration contours of OH species and H 2 0 
species, respectively. The contours indicate reacted regions as expected. 
The line plot in Fig. 4.4 illustrates distributions of various species along 
the y-location of about 0.13 cm from the lower wall (y-station No. 15). It 
is seen that the OH concentration increases rapidly across the shock wave 
and stays relatively constant thereafter (at its equilibrium value). The 
H 2 0 concentration, however, increases continuously signifying that more fuel 
is being consulted by an exothermic process. The behavior of the species 
distributions is in good accordance with the physics of the global two-step 
chemistry model used. A similar plot at the lower-wall surface is given in 
Fig. 4.5. 

The geometry of a more realistic problem and its inflow conditions are 
shown in Fig. 4.6. The geometry and flow conditions are representative of a 
planar cut perpendicular to the inward swept leading edge of a three- 
dimensional single strut scramjet engine. It is noteworthy that this is the 
same configuration used in the study of Drunmond and Weidner [74]. Gaseous 
hydrogen fuel is injected into the flow from the side walls and from tne 
center strut such that the overall equivalence ratio is equal to one. The 
conditions of injected hydrogen are: p = 258706 N/m 2 , T = 246 K, u = 156 
m/sec, v = 1241 m/sec. The injectors are located 6 cm downstream of the 
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Fig. 4.5 Distributions of species mass fraction at the lower wall. 
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Geometry and flow conditions of case No 



engine cross sectional minimum and are each 0.1 cm wide. Due to symmetry, 
only the upper half of the domain is solved. A grid system of 87x30 is 
employed. Grid points are concentrated near the solid surfaces in tne y- 
direction and near the injectors in the x-direction (see Fig. 3.2). Clus- 
tering of mesh points in the near field of the injectors is necessary to 
properly resolve jet-flow interaction and the mixing process between the 
fuel and air. The computational grid used is the same as that of Ref. 74 in 
order to allow point by point comparisons of the results. 

It is of interest to study the nonreacting (mixing only) case first so 
that the inferences can be made later on as to the effects of combustion on 
the flow field. Figures 4. 7-4.9 show the pressure contours, velocity 

vectors, and H 2 concentration contours, respectively, for the nonreacting 
case. The figures reveal that small separation bubbles occur before the 
throat due to interactions of shock waves and boundary layers. Separations 
also occur immediately upstream of the injectors; it is imperative to have 
these separations for proper flame holding in the engine. 

Results with chemical reactions are given in Figs. 4.10-4.15. Figures 
4.10-4.12 illustrate the similar information to that in Fig. 4. 7-4. 9. It is 
seen from Fig. 4.11 that the separation bubbles upstream of the throat 
increases in size as compared to the nonreacting case. This is believed to 
be due to a higher level of back pressure caused by the heat release from 
combustion. The separation bubbles, in turn, induce upstream shock waves as 
can be seen in Fig. 4.10. The occurrence of these bubbles due to inlet- 
combustor interaction can easily choke the flow in the inlet and should be 
avoided, if possible, in design. The hydrogen contours shown in Fig. 4.12 
indicate that the fuel jets do not penetrate significantly into the core air 
flow resulting in low level of mixing. Increasing the jet pressure and 
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Fig. 4.15 Upper surface pressure distributions. 
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velocity can, again, choke the flow. Parametric studies of the jet condi- 
tions required to promote more mixing will be the subject of future study. 
Figure 4.13 shows the water concentration contours of the reacting case 
indicating the regions where combustion has taken place. 

One major purpose of this study has been to compare the extent of com- 
bustion using the present finite rate chemistry model with the results of 
Drummond and Wiedner [74] where a complete reaction model was used. In the 
complete reaction model, only the species equation for the fuel is solved 
along with the fluid dynamic equations. There is no source term in the fuel 
species equation and therefore the system is not stiff. The extent of com- 
bustion is determined by assuming that when the fuel and oxidizer mix at 
temperatures above 1000 K they react to their maximun extent possible 
(stoichiometric limit) until one reacting partner is completely depleted. 

The comparisons of water profiles at some selected x- stations are made 
in Fig. 4.14. The cross stream coordinate at each station was normalized 
with the corresponding channel's width such that its value ranges from zero 
to one. The top graph shows the profiles upstream of the injectors within 
the recirculating zones where flames are stabilized. It is seen that the 
concentration levels of water, which are an indication of the extent of 
combustion, are approximately the same. This could be attributed to the 
fact that mixture velocity in the recirculating zones is relatively small 
providing ample resident time for it to react to an appreciable amount. The 
bottom graph of Fig. 4.14 indicates the water profiles downstream of the 
injectors where it can be seen that the finite rate model predicts a lower 

level of combustion than the complete reaction model. This is consistent 
with the formulation of the two models since the finite rate schane does not 
necessarily allow complete reaction. 
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Finally, normalized surface pressure, (P-P max )/(P max -P m i n ) » is plotted 
in Fig. 4.15. The highest spike in the plot indicates the location of the 
injector. The finite rate model predicts a small bump in surface pressure 
upstream of the cross sectional minimum whereas the complete reaction model 
does not. This pressure bump was caused by the induced shock wave ahead of 
the separation bubble; the complete reaction model does not show as large a 
separation. Downstream of the injector, however, the pressures exhibit 
similar trends except that the pressure predicted by the complete model is 
somewhat higher due mainly to a greater degree of combustion. 

4.2 Results of. Parabol ized Code 

All results reported in this section were obtained by using the single 
sweep method with the value of u> equals to 0.7. In most cases, the values 
of important parameters used are indicated at the bottom of the figures. 

The premixed test case (Fig. 4.1) solved earlier by the elliptic code 
has been reinvestigated by the parabolized code. Initial data profiles were 
obtained from those of x-station No. 5 of the elliptic results which is 
located at five stations upstream of the compression corner. Regardless of 
the schemes used, the code is always started up by the Euler implicit scheme 
which requires only one level of initial data. In the following discussion, 
comparisons are made between results of the parabolized and elliptic codes. 

Results obtained by using the second order accurate three point back- 
ward scheme are illustrated in Figs. 4.16-4.18. Comparisons of temperature 
and pressure between the parabolized (PNS) and elliptic (NS) results are 
made in Fig. 4.16. The plots are made along the y-station located approxi- 
mately 0.13 cm from the lower wall (y-station No. 15). It can be seen that 
the comparisons are very reasonable. In the absence of combustion, the 
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4.18 Cross stream profiles at the outflow boundary fo 
various variables, three-poinc-backward scheme. 



pressure and temperature behind the shock wave would remain constant. 
Increase in pressure and temperature due to heat release during combustion 
are predicted very well by the PNS code. Both NS and PNS results show some 
oscillations in pressure. Use of parameter C c = 0.2 and/or larger values of 
the smoothing coefficients do not suppress these oscillations. The distri- 
butions of H 2 0 and OH species along the same y-station are compared next in 

Fig. 4.17. The H 2 O mass fractions compare well for both formulations but 
the OH mass fractions show some disagreement. In fact, OH mass fraction in 
the PNS results at the grid point right at the shock front takes a huge jump 
out of the plot frame. Fortunately, the jump occurs at only one grid point 
and appears to be a very local phenomenon. By refining the mesh spacing in 
the streamwise direction two times, however, the jump disappears as indi- 
cated in Fig. 4.19. It is, therefore, concluded that the jump in OH species 
(which is the one that causes chemical stiffness) occurs because of large 
linearization errors due to the use of a coarse mesh. Even with the refined 
mesh, the PNS results underpredict the OH mass fraction when compared with 
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lagged temperature used in the evaluation of the Jacobian matrix for the 
chemistry source terms. As the flow develops further downstream, however, 
OH concentration picks up to the same level as the NS results. Figure 4.18 
details the comparisons of various quantities across the channel at the 
outflow boundary of the geometry. The cross stream coordinate, Y, was nor- 
malized such that its value ranges from zero to one whereas the dependent 
variables are scaled with respect to their corresponding maximum values. 

All profiles appear to be in reasonable agreement. Both results suffer 
oscillation in pressure across the shock front which is typial of a central 
difference scheme. Not shown here are the velocity profiles; numerical 
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output indicates that velocity profiles for the two formulations are in good 
agreement. 

Figure 4.20-4.22 show the same information as those of Figs. 4.16-4.18 
but the new figures were produced by using the first order accurate Euler 
implicit scheme. Note that the value of C c is equal to 0.2 rather than 
zero. The vanishing of C c for this scheme cannot be allowed since it gives 
rise to an instability regardless of the magnitude of the values of the 
smoothing coefficients; this appears to be the case only for reacting flows 
since zero value of C c was used all the time in nonreacting flow calcula- 
tions [8, 9]. A jump in OH concentration across the shock front persists as 
in the three point backward results. In general, these plots indicate that 
results of the Euler implicit scheme are slightly inferior to those of the 
three point backward scheme. 

The single strut scramjet engine problem is also reinvestigated with 
the PNS code. The x-station No. 59 was selected as the starting solution; 
this is reasonably well downstream of the injectors where elliptic effects 
should be weak. Like the previous case, only the one-level Euler implicit 
scheme is used to start up the integration. In all the plots to be discuss- 
ed, the symbol always indicates the PNS results whereas the solid line 
represents the NS results. Figures (a) always indicate surface pressure . 
whereas Figs, (b), (c) and (d) always indicate cross stream pressure at the 
outflow boundary (x-station No. 85), cross stream pressure at x-station No. 
70 and cross stream mass fractions of H 2 O at x-station No. 70, respective- 
ly. All variables are scaled with respect to their corresponding maximum 
values and the spatial coordinates are normalized such that their values 
range from zero to one. 

Results of the three point backward scheme are shown in Fig. 4.23. It 
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Fig. 4.22 Gross stream profiles at the outflow boundary for 
various variables, Euler implicit scheme. 


80 






is seen that the early phase of surface pressure (Fig. 4.23(a)) is oscilla- 
tory. The cross stream variables at x-station No. 70 compare very well with 
the NS profiles but the cross stream pressure at x-station No. 85 is not in 
good agreement. By using C c = 0.2, the oscillation in surface pressure 
disappears as can be seen in Fig. 4.24(a). However, quality of other vari- 
ables is slightly degraded; this is probably due to the fact that a non- zero 
C disrupts the second order accuracy of the scheme. Larger values of C 
worsen the PNS solution even further but the surface pressure becomes 
smoother. Increasing the values of C. and C F , while using C = 0, however, 
appears to improve solution of the PNS code as is evident from Figs. 4.25 
and 4.26. 

Results obtained by using the Euler implicit scheme are presented in 
Figs. 4.27-4.29. As in the previous problem,' the scheme becomes unstable 
for the value of C c equals to zero. In general, the results shown in these 
three figures are superior to the results of the three point backward 
scheme. This surprising paradox was reported earlier by Chitsomboon et al . 
[S] is solving nonreactiong flows. It appears that the Euler scheme has 
more natural damping capability embedded in it since the first truncated 
term in the scheme is an even order term (second-order derivative) whereas 
the three point backward scheme has an odd order term (third-order deriva- 
tive) as its first truncated term. An even order derivative term, with 
proper sign, behaves like an artificial viscosity terms which helps in damp- 
ing disturbances whereas an odd order derivative term disperses the disturb- 
ances without damping their magnitudes [38]. 

Figure 4.27 shows the results by using C c = 0.2. The cross stream 
pressure at x-station No. 70 and the surface pressure are slightly dissi- 
pative in comparison with their three point backward counterpart in Fig. 
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Fig. 4.25 Results of three-point-backward scheme (No. 3). 
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4.23. The pressure at x-station No. 85, however, compares more favorably 
with the NS results. The solution is further improved by increasing C c to 
0.5 as is seen in Fig. 4.28. Larger values of Cj and do not seen to 
affect the quality of the solution (Fig. 4.29). 

5. CONCLUDING REMARKS 

Two computer programs for solving viscous supersonic chemically react- 
ing flow problems associated with the hydrogen-air system have been develop- 
ed. In the first computer code, the unsteady form of the fully elliptic 
governing equations are anployed in order to capture strong upstream inter- 
actions of the flow. The finite rate global two-step chemistry model is 
used to represent hydrogen-air combustion. The stiffness resulting from 
using the chemistry model is circumvented by evaluating the chemistry source 
terms implicitly. The elliptic code is written in a vector FORTRAN language 
which is made operational only on the VPS-32 vector processing computer at 
the NASA-LaRC. The fluid dynamics part of the governing equations is 
integrated in time by a fully explicit MacCormack algorithm. Due to implic- 
it source terms, the chemistry part of the governing equations gives rise to 
a block monodiagonal system of algebraic equations which can be solved also 
in a vectorized manner giving considerable saving in computer time over a 
scalar inversion procedure. The results obtained from the elliptic code in 
solving two internal flow problems are physically reasonable. In the case 
of the scramjet inlet-combustor configuration problem, the code predicts 
strong interactions between the inlet and the combustor which could choke 
the flow in the inlet. 

The second computer code employs a parabolized form of the governing 
equations in which only a fraction of the streamwise pressure gradient is 
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retained. The code is applicable for supersonic high Reynolds number flows 
with weak upstream interaction such as regions far downstream of the trans- 
verse fuel injector in a scramjet engine. The governing equations are solv- 
ed in a fully coupled fully implicit manner in order to capture strong shock 
waves as well as fluid-chemistry interactions. The finite difference algo- 
rithm used is fairly general; it degenerates into many familiar schemes upon 
specifying two parameters. The fully coupled procedure yields a Dlock (8x8) 
tridiagonal system of algebraic equations which is solved by a noniterative 
method. Solution is obtained by forward marching in the streamwise spatial 
direction only one time providing considerable saving in computer resources 
when compared with the elliptic code. If desired, a multi-sweep procedure 
can be devised with minimal effort. The results obtained from the para- 


hoi ized code cofnpdre redsondbly well with those from the elliptic results* 
The first order accurate Euler implicit scheme is found to give good 
results. In some cases, the Euler implicit scheme gives even better results 
than the second order accurate three point backward scheme. In a real 
application, it is highly desirable to have a bench mark solution to compare 
results of the three point backward scheme against results of the Euler 
implicit scheme in order to determine which scheme is more appropriate for 
the class of flow being investigated. For cases where bench mark solutions 
are not available, the Euler implicit scheme is recommended since, in 
general, it should provide fairly accurate results. If the three point 
backward scheme is selected, the suggested values of various parameters are: 
Cj = 0.25, = 0.125, C c = 0.0. For the Euler implicit scheme, the follow- 
ing values are recommended: C T = 0.25, C c = 0.125 and C = 0.5. 

i t c 

The robustness of the parabolized code is considerably improved by 
introducing the parameter C c into the continuity equation and the species 
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equations. The value of C * 0.2 can suppress pressure oscillation while 

V 

still maintaining overall accuracy. Non-zero value of C c , however, is 
recommended for a first order accurate scheme only since it could disrupt 
the formal accuracy of higher order schemes. 

The results of the two problems solved by the present parabolized algo- 
rithm confirm that the algorithm is highly conservative. Regardless of the 
quality of the solution, mass flow rate across the channel at all stations 
are always within 0.01 percent of one another. 

Further investigations are obviously needed in the area of turbulent 
jet mixing. The algebraic turbulence model used in this study serves only 
as a starting point. The multi-component mass diffusion processes also need 
more elaboration. In any event, experimental data are highly desirable to 
validate the models. More sophisticated finite rate chemistry models can be 
included especially in the parabolized formulation. Inclusion of better 
chemistry models always means more species and more reaction paths; this 
makes the elliptic equation approach very expensive. 
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APPENDIX A 


DERIVATION OF HEAT FLUX TERMS 

Only the 3q/3x term will be derived here; relation for 3q/3y follows 
the same line of derivation. The heat flux term is given in its unsimpli- 
fied form as 


. _ 5 3f . 

q - - k II - pD E [(— L) h ] (A.l) 

3x i=l 3x 

by using the relation 

a - k/pCp (A. 2) 

and 

Le 3 a/D . (A. 3) 


Eq. (A.l) can be rearranged as 

q 3 -pD{LeC — + Z [(— ) h .]} (A.4) 

* p 3x i=l 3x 1 

with the help of the relations (A. 3) and 

Pr 3 v/o (A. 5) 


the term pD can be written as 
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pD = 


Le Pr 


For a unit Lewis number, Eq. (A. 4) can be given as 


q x - - — [c — + i {(i_ f )n J] 

Le Pr P 3x i=l 3x 

where 


C = z f. C • 
P i=i 1 P 1 


The static enthalpy for the mixture is given by 


= j r i f, («; * / c p( u) 




Note that the dummy variable, 5, is employed in evaluating the 
enthalpy. Upon taking a partial derivative with respect to x < 
the Leibnitz rule one obtains 


3h 


3f . 3H°. 


(( H l + L 


3x i=l 


3x 3x 


T 3C . 

f _ F 

0 3x 


+ f / —El (0 d? + f C (T) II] 

I r\ * „ 1 P I <\ 


3X 


By realizing that the coefficient of the first term is nothing 


(A. 6) 


(A. 7) 
(A. 8) 


(A. 9) 

sensible 
,nd applying 


(A. 10) 
but h . and 
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the second and the third terms are identically zero, Eq. (A. 10) can be 
written as 


9h 

9x 


5 

E 

i = l 


9f. 5 

( ) h. + e 

9x 1 i=l 


f. c .11 

1 P 1 9x 


(A. 11) 


or 


ah 5 9f. T 

— = £ ( — L) h, ♦ r. — • (A. 12) 

9X i=l 9x w 9x 

Substituting Eq. (A. 12) into (A. 7), the desired relation is finally obtain- 
ed. It is important to note that the simplification is possible only for 
the unit Lewis number assumption. The simplification is useful since it 
helps in organizing and simplifying the coding chore and it offers some 
saving in computer storage because only one array, h, is needed instead of 
five arrays of h^ . 
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APPENDIX B 


JACOBIAN MATRIX FOR CHEMISTRY SOURCE TERMS 
Components of the Jacobian matrix for the chemistry source terms, M.., 

J 

are given in this appendix. The suffixes i and j denote row and column of 
the matrix, respectively. Whenever a component is not given, it means that 
component is identically zero. It should be noted that temperature depen- 
dency of the kinetic rate terms, and k b , is not considered in evaluating 
the components. The temperature for these terms is simply lagged one step 
behind the step where solutions are being sought. The matrix takes exactly 
the same form for both the elliptic and the parabolized formulation out in 
the parabolized formulation the matrix is defined as C rather than M (see 
section 3.2). The components are given as follows: 

z - J n+1 /J n 


Mil = -(k fl /M 3 )Zp 3 
Ml 3 = -(k fl /M 3 )Z p x 
M14 = (k bl M^MjjZ p 4 


M22 = -(4k b2 /M 2 )Z p 2 


M2 3 


2k f2 M 2 Z 
M 3 H 4 


2 

p 4 
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4k f , M, Z 
M2 4 = L£ £ 

M 3 M 4 

Cl = M /M 
C2 = -M 3 /2M 2 
M31 = (Cl) (Mil) 
M32 = (C2) (M22) 
M33 = (Cl) (M13) • 
M34 = (Cl) (M14) 
C3 = -2M./M, 

H i 

C4 = -M 4 /rt 2 
M41 = (C3) (Mil) 
M42 = (C4) (M12) 
M43 = (C3) (M13) 
M44 = (C3) (M14) 


(C2) (M23) 
(C2) (M24) 


(C4) (M23) 
(C4) (M24) 
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APPENDIX C 


VECTORIZED SUBROUTINE FOR SOLVING 
A BLOCK MONO-DIAGONAL SYSTEM 

The following subroutine is written in CDC CYBER 200 FORTRAN VERSION 2 
(Vector FORTRAN) in order to make it operates efficiently on the VPS-32 
vector processing comuter at the NASA-LaRC. The solution procedure can be 
made more efficient by realizing that only a fraction of mesh points are 
chemically reacting; the rest of the mesh points are nonreacting and can be 
integrated by a fully explicit procedure which is much less expensive. The 
reacting mesh points are distributed throughout the domain, more or less, 
randomly. The scalar array, (IT), and the bit array, (B2), are calculated 
independently in the main body of the program and are used as maps for 
determining which mesh points are chemically active. The sparse arrays can 
then be compressed into concatenate locations giving shorter vectors lengths 
and therefore less mathematical operations. Finally, the solution can be 
expanded back into their real locations with the help of the maps (IT) and 
(B2) . 


SUBROUTINE LUDEC(NN) 

BIT B1,B2 

C0MM0N/C0M1/IT(NY,NX),B2(NY,NX),Q(NY,NX,4,4),RHS(NY,NX,4) 
1,SLN(NY,NX,4) ,EBAR(NY,NX,4) ,NXM4, IRP 
C 

C DECOMPOSE THE ORIGINAL (Q) MATRIX TO (LU) VECTORIALLY. 

C THE METHOD IS A STRAIGHT FORWARD DECOMPOSITION WITH THE 
C 0IAG0NAL TERMS OF (L) EQUAL TO ONE. 

C THE (LU) MATRIX IS STORED BACK INTO THE (Q) MATRIX. 
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C NOTE ALSO THAT THE SPARSE VECTORS ARE COMPRESSED INTO 
C CONCATENATE VECTORS, BY THE FUNCTION "Q8VCMPRS\ FOR 
C EFFICIENT OPERATION. THE COMPRESSION AND EXPANSION BACK 
C INTO THE SPARSE VECTORS ARE CONTROLLED BY THE ARRAY 

C (IT) AND THE BIT VECTOR (B2) 

C (RHS) IS THE UNKNOWN 

C (SLN) IS TEMPORARY ARRAY 

C (Q) IS THE ORIGINAL BLOCK MONO-DIAGONAL MATRIX 
C IRP IS THE LENGTH OF CONCATENATE VECTORS 

C NN IS THE DIMENSION OF THE BLOCK MATRIX 

C 

NXM4=NX*NY-2*NY-2 
DO 5 1=1, NN 

5 WHERE (IT(2,2; NXM4 ) . EQ.O) SLN (2, 2, 1 ;NXM4)=RHS(2,2, 1 ;NXM4) 
C 

C COMPRESS THE RHS 
C 

DO 6 1=1,4 

RHS(2 ,2, I ;NXM4)=Q8VCMPRS(RHS(2 ,2, I ;NXM4) ,B2(2 ,2 ;NXM4) ; 
&RHS(2,2, I ;NXM4) ) 

6 CONTINUE 

DO 100 1=2, NN 
DO 110 J=1,NN 
IF( J-I )10,20,20 
10 M=J-1 

IF(M)30,30,40 
40 DO 50 K=1,M 
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50 Q(2,2,I,J;IRP)=Q(2,2,I,J;IRP)-Q(2,2,I,K;IRP) 

&*Q(2,2,K, J; IRP) 

30 Q(2,2,I , J;IRP)=q(2,2,I, J;IRP)/Q(2,2,J, J;IRP) 

GO TO 110 
20 M= I - 1 

DO 60 K=1,M 

60 Q(2,2,I,J;IRP)=Q(2,2,I,J;IRP)-Q(2,2,I,K;IRP) 

$*Q(2,2,K, J; IRP) 

110 CONTINUE 
100 CONTINUE 
C 

C SOLVE (LU)X=RHS 

C X IS STORED BACK INTO RHS AND LATER ON TO SLN 
C THE LAST STEP IS NOT NECESSARY, I.E., THE FINAL RHS 
C CAN BE USED AS THE FINAL SOLUTION. 

C 

C SOLVE (L)X*=RHS 
C 

105 DO 200 M=2,NN 
N=M-1 

DO 140 K=1,N 

140 RHS(2,2,M;IRP)=RHS(2,2,M;IRP)-RHS(2,2,K;IRP)*Q(2,2,M,K;IRP) 
120 CONTINUE 
200 CONTINUE 
C 

C SOLVE (U)X=X* 

C 
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DO 300 M=NN,1,-1 
IF(M-NN)320,310,310 
320 N=M+1 

DO 330 K=NN,N-1 

330 RHS(2,2,M;IRP)=RHS(2,2,M;IRP)-RHS(2,2,K;IRP)*Q(2,2,M,K;IRP) 

310 RHS(2,2,M; IRP)=RHS(2,2,M; IRP)/Q(2,2,M,M; IRP) 

300 CONTINUE 
C 

C EXPAND BACK THE SOLUTION INTO (SLN) 

C 

DO 340 1=1,4 

EBAR(2,2, I ;NXM4)=Q8VXPND(RHS(2,2, I ;NXM4) ,82(2,2 ;NXM4) ; 
&EBAR(2,2, I ;NXM4) 

WHERE ( IT(2,2; NXM4) . EQ. 1 ) SLN (2,2, I ;NXM4)=EBAR(2 ,2, I ;NXM4) 
350 CONTINUE 
RETURN 


END 


